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ABSTRACT 



How to create planetesimals from tiny dust particles in a proto-planetary disk 
before the dust particles spiral to the central star is one of the most challeng- 
i ng problerns in th e theory of planetary system formation. In our previous paper 
kato et ahlbopgh . we have shown that a steady angular velocity profile that con- 



sists of both super and sub-Keplerian regions is created in the disk through non- 
uniform excitation of Magneto- Rotational Instability (MRI). Such non- uniform 
MRI excitation is reasonably expected in a part of disks with relatively low ion- 
ization degree. In this paper, we show through three-dimensional resistive MHD 
simulations with test particles that this radial structure of the angular veloc- 
ity indeed leads to prevention of spiral-in of dust particles and furthermore to 
their accumulation at the boundary of super-Keplerian and sub-Keplerian re- 
gions. Treating dust particles as test particles, their motions under the influence 
of the non-uniform MRI through gas drag are simulated. In the most favorable 
cases (meter-size dust particles in the disk region with a relatively large fraction 
of MRI-stable region), we found that the dust concentration is peaked around 
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the super/sub-Keplerian flow boundary and the peak dust density is 10,000 times 
as high as the initial value. The peak density is high enough for the subsequent 
gravitational instabihty to set in, suggesting a possible route to planetesimal 
formation via non-uniformly excited MRI in weakly ionized regions of a disk. 

Subject headings: protoplanetary disks — instabilities — MHD — planetary 
systems: formation — turbulence 



1. Introduction 



One of the most serious problems in planet formation is how the dust components 
grow up to larger objec ts in protoplanetary disks. The kn own key-stage is when the dust 
particles are meter-size (lAdachil Il976l : IWeidenschilling] 119771 ) . Since gas pressure gradient in 
a protoplanetary disk is negative, gas rotates at sub-Keplerian velocity. On the other hand, 
the dust particles tend to rotate at Keplerian velocity and therefore they feel headwind of 
gas. Losing their angular momenta, they spiral into a central star quickly on timescales of 
100 years. 

One possibility to overcome the meter-size barrier is the pla netesimal formation via grav- 

itational collapse of a cluster composed with small dust particles dSafronovlEoGoi iGoldreich fc WardI 
19731 ). The dust would settle to midplane in the absence of turbu lence. Earl y stud - 



ies of dust settling i n turb ulen t accretion disks wer e carried out by ICuzzi et al.l ( 1l993l ). 
Miyake fc Nakagawal (119951 ) and iDubruUe et al.l ( 119951 ) . They presented that large dust par- 
ticles settle down toward their equilibrium distribution in a turbulent diffusive time scale 
while small ones remain mixed throughout the whole gas disk. If dust density in midplane 
becomes high enough, Kelvin-Helmholtz instabilities between a sediment layer of dust and 
the overlying gas c an drive turbulence and mi x the dust layer enough to prevent gravi- 



tation al instabilitv (IWeidenschilhnd Il977l Il980l : llshitsu fc Sekival 12003 



Gomez fc Ostriker 



20051 ). iJohansen et al.l (l2006bl ). however, carried out two-dimensional simulation and found 
active turbulent concentration. Except the effect of Kelvin-H elmholtz instabilities, t he 



high dust density iii midp lane leads to a streaming instability (iGoodman fc Pindorl 12000 



Youdin fc Goodmanll2005l) because the dust is susceptible to preferential radial migration rel 
ative to the gas (INakagawa et al.lll986l ). The streaming instability was shown to lead to tur- 
bulen ce and concentrate dust particles locally ( lYoudin fc Johansenll2007l : IJohansen fc Youdin 
20071 ). As one of other mechanisms of turbulent con centration, the effect of t he turbulence 



drive n by the magnetorotational instability (MRI; iBalbus &: Hawleyl Il991t iHawley et al. 



19951 ) is recognized. P lanetesimal formation by gra vitati onal instability in t he MRI tur- 
bulence is addressed by IJohansen et al.l (l2006al . 120071 ) and iBalsara et al.l (120091 ) . 
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Turbulence level is the most fundamental parameter governing whether planetesimals 
can form by gravitational instability in these models. The turbulence level of MRI depends 
on the ionization degree of disk gas and the strength of magnetic field. The linear analysis 
showed that the growth rate of MRI is reduced by ohmic diffusion when collisions are so 
frequent that not only the ions and neutrals are well coupled but also electrical currents 
are damped (jjinlll996l : ISano fc Miyamal 119991 ). At the same time, the growth wavelength 
becomes longer and transport rate of angular momentum is diminished. The nonlinear 
evolution of the MRI in non i deal MHD starting from the weak vertical magnetic field has 
been studied by ISano et al.l ( 119981 ). They found that sustained MHD turbulence requires 
the magnetic Reynolds number = v\/riQ > 1.0, where va, V and Q are Alfven velocity, 
magnetic resistivity and Keplerian rotation frequency. 

The values of rj are regulated by an ioni zation structure in a protoplanetary disk. Th e 



inner regions of a disk are thermally ionized (jPneuman fc Mitchell 



19651 : lUmebavashil 119831 ). 



At greater distances from the disk center, stellar X -rays and diffuse cosmic rays ionize the 
surface gas layer down to a certain column density (iGlassgold et al.l 119971 : llgea fc Glassgold 
19991 ). In moderately distant regions in which disk column density is high enough, midplane 
layers are not highly ionized and MRI could be inactive there. This region is called "dead 
zone" (jGammidll996l ). 



The characteristics and interesting physics are prospective around dead zones. The 
magnetic turbulence outside the dead zone, which transports angular momentum and mass, 
creates gas density bump at the outer boundary of the dead zo ne. This density pile-up 
triggers the Ro ssbv wave (iLi et al. IVarniere fc Tagged 120061 ) to form long-life anticy- 

clonic vortices. iBarge fc Sommeria (1995) suggested that the dust particles are trapped in 
a center of vor t ex. B y numerical simulation of Rossby vortices includi ng solid pa r ticles , 
Inaba fc Bargd (120061 ) confirmed the enhancement of dust density there. iLyra et al.l (120081 ) 
showed that gravitationally bound embryos can form in the regions of enhanced dust density. 



Kato et al.l (|2009l . ; hereafter referred to as Paper I) showed that non-uniform excita- 



tion of MRI creates angular velocity profile that consists of super and sub-Keplerian regions 
locally in the disk (for details, see section 2). It is expected that this radial structure of 
the angular velocity leads to prevention of spiral-in of dust particles and their accumulation 
at the boundary of super and sub-Keplerian regions, since the super-Keplerian flow pushes 
the dust particles outward while the sub-Keplerian flow drags them inward. This interesting 
mechanism is based on the co-existence of MRI active and inactive regions. The ioniza- 
tion de gree that regulates MRI depends on t he relative abundance and size distribution of 
grains (jWardle fc Ng Il999l : ISano et al.l l2000l ). It can also be changed by phase chemistry 
(lllgner &: Nelson! |2008| ). Slight changes in these quantities near boundaries at global ac- 
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tive/dead zones of MRI may create the local inhomogenity in MRI growth. We will address 
this issue in a subsequent paper. Note that the local model in this paper is not a toy model 
that mimics the global inhomogeneous structure of dead and active zones. 

In this paper, we investigate dust accumulation induced by inhomogeneous MRI evo- 
lution by three-dimensional non-ideal MHD simulation including dust particles. The dust 
particles are represented by Lagrangian particles suffering gas drag and are assumed not to 
affect on gas. We consider both cases with radially non-uniform magnetic resistivity under 
uniform magnetic field direction and non-uniform magnetic field with radially uniform mag- 
netic resistivity. MRI develops in low resistivity regions in the former case and in strong 
vertical magnetic field regions in the latter case. The latter model is the same as the model 
adopted in Paper I. We briefly summarize the results of Paper I in section [21 We explain our 
simulation setting in section [31 In section [H we present that the MRI leads to quasi-steady 
state and particles are locally concentrated as expected. We also analyze the distributions of 
particle radial velocity to understand the effect of weak remnant magnetic turbulence. We 
examine how high the particle density enhancement is and the possibility of planetesimal 
formation via gravitational instability. Finally, we summarize our results and discuss their 
implications in section [SI 



2. Brief Summary of [Paper I 



In Paper I, we have studied how the angular velocity profile of gas is modified when 
MRI is excited non-uniformly in a part of a disk. Local two-dimensional simulations (in the 
x-z plane where x and z are radial and vertical coordinates to the disk) with the shearing 
box boundary conditions were performed. We assumed the ionization degree of the gas to 
be relatively low. The vertical component to the disk plane {z component) of magnetic field 
was changed along the radial direction (x) while the absolute values of the field is constant. 
In the part of a disk where the ionization degree is low (equivalently, resistivity is high), MRI 
is active only in the regions with relatively large vertical component of the magnetic field. 
One of the models used in the present paper has the same setting but with three-dimensional 
simulations, while the other model in the present paper assumes non-uniform resistivity with 
uniform magnetic field. 

The most illustrative result of Paper I is shown in Figure [H Time evolution = 
0,40,70) of vertically averaged pressure and angular velocity of gas is plotted in Panels a 
and b. MRI is excited only in the initially unstable region (—0.71 < x/H < 0.71) with large 
enough initial vertical field and vigorous angular momentum and mass transport occurs 
within the unstable zone. When the magnetic perturbations are propagated to the initially 
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stable regions {\x/H\ < 0.71), they are dissipated via resistivity and the large part of the 
stable region remains undisturbed. Hereafter, we call the initially MRI stable and unstable 
regions simply as "stable" and "unstable" regions, respectively. After the vigorous angular 
momentum and mass transfer, a new quasi-steady state with the local rigid rotation feature 
appears in the unstable region (Figure [TJ3). The pressure profile of the gas is also modulated 
(Figured^). The rigid rotation feature that appears in the saturated quasi-steady state is 
sustained by the modified pressure pattern with its two peaks located around the unstable- 
stable boundary regions and its bottom located in the middle of the unstable zone. The 
quasi-steady rigid rotation pattern results in super-Keplerian angular velocity at the outer 
part of the unstable region. Since the gas in the stable region is unperturbed and remains 
sub-Keplerian when the global pressure gradient is taken into account, one finds that super- 
and sub-Keplerian zones co-exist in the part of the disk due to the non-uniform excitation 



The co-existence of super- and sub-Keplerian zones has significant implication for the 
dust particle dynamics in a disk. As mentioned in the Introduction, the dust particles feel 
headwind and migrate inward in the sub-Keplerian gas flow. On the other hand, the dust 
particles feel tailwind and migrate outward in the super-Keplerian zone. Then at a boundary 
between sub- and super-Keplerian zones, with the super-Keplerian zone situated on the inner 
side, which is indeed seen at the outer edge of the unstable region, the inward migration of 
dust particles through the sub-Keplerian zone (the stable zone) are halted. The boundary 
region potentially acts as a barrier for the dust infall and furthermore as the focal annulus 
of the dust particle concentration which may eventually lead to planetesimal formation. 

The widths of the stable and unstable regions are one of the key parameters that deter- 
mine the characteristics of the final state. We have found that the final state is classified by 
the spatially averaged magnetic Reynolds number in the initial condition. 



where v^z = B^/ y/^T^p is z component of Alfven velocity. If -Rm.ave is less than unity, a 
quasi-steady state with the local super-Keplerian rotation appears. The turbulence level 
within the stable region and at the boundary varies according to -Rm.ave- The degree of 
dust concentration at the boundary region would be affected by the level of turbulence. For 
-Rm,ave > 1, the magnetic perturbations are not dissipated enough within the stable region 
and the velocity pattern undergoes perpetual changes without reaching a quasi-steady state. 
Steady accumulation of dust particles is not expected in this case. Note that our simulations 
with constant -Rm,ave converge with increasing resolution, because the resistivity 77 (ohmic 
dissipation) is explicitly included in our simulations and it is kept constant due to the 

constancy of -Rm.ave- 



of MRI. 



■m,ave 



^A^.ave/^^, 




In the present paper, we address the following questions that arise from the results of 
Paper I: (1) Do dust particles truly concentrate at the boundary? (2) Is the emergence of 
the quasi-steady rigid rotation region preserved in three-dimensional cases? (3) How does 
remaining weak turbulence affect the dust concentration process? (4) Are dust particles 
concentrated dense enough on long enough timescales for subsequent gravitational instability 
to set in? 



3. Model 



3.1. Equations 



Here we are concerned with a small region in a protoplanetary disk and study the 
local dynamics of gas and dust particles. The gas motion is calculated by resistive MHD 
simulations. The dust particles are modeled as test particles that move under the effects of 
the gravity from the central star and the gas drag. 

We adopt a local shearing box that has radial (x), azimuthal (y) and vertical (z) di- 
rections. The boundary conditions are periodic in the y and z directions and the shearing 



box b oundary condition is used in the x direction (jWisdom &: Tremainelll988t iHawley et al. 
19951 ). We calculate compressible resistive magnetohydrodynamic (MHD) equations given 
by 



- + (u.V)u 

dB 

'dt 
P 



P \ Stt 

0, 



47rp 



;B ■ V) B - 212 X u + Sfi^xx - /3c,^]x, (2) 



V X [(u X B) -r/(V X B)] , 



cIp, 



(3) 

(4) 
(5) 



where x is a unit vector in the x direction. We consider only the Ohmic dissipation in equa- 
tion (jlj) because influence of ambipolar diffusion is negligible at t he density level that is ex- 
pected in the inner (< lOOAU) mi d p lane of a protoplanetary disk ( 1jinlll996l : ISano fc Miyama 
19991 : IChiang fc Murray-Clayll2007l ). The last terms in r.h.s. of equation ([2]) allows the effect 
of the global pressure gradient to be taken into account in our local model. The parameter 
(3 is related to the global pressure gradient as 
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aCs^ 

r 



(7) 



\l r 



a 



where we assume P oc r", with constant a and an isothermal disk = HQ, where H is disk 
scale height. The value of P (= aH/r = acg/rfl) is treated as a constant paramet er in our 
local model. We solve equation @ with MOCCT method (IStone fc Normanlll992l ) and the 
advection terms in equation ([2]) with CIP method (lYabe fc Aokilll99ll ). 



The dust grains are treated as test particles. The equation of motion of the i-th particle 
is given by 



^ = -2r2 X Vj + 31]^XjX - — (vi - Wi) , 
dt Tf 



(9) 



where Wj is the gas velocity at the location of the z-th particle, which is interpolated with 
the three-dimensional first-order interpolation scheme by values of u at the eight grid points 
surrounding the particle. 

The friction force, the last term in equation (j9]), is proportional to the velocity difference 
between the dust and gas since the dust sizes that we are concerned with are in Epstein and 
Stokes drag-law regimes. The characteristic friction (stopping) timescale tj is given by 



alp. 



min {a,Cs, 9/2z/) p' 



(10) 



where a,, p, and Cs are the dust particle radius, dust particle density and the sound speed, 
respectively, u is the molecular viscosity and equal to CsA/2, where A is the mean free 
path of the gas molecules. Depending on the particle size relative to the mean free path 
(a, > or < (9/4) A) the drag law changes from the Epstein to the Stokes regime. The scaled 
friction times of TfQ = 0.1 and 1.0 correspond to the grain sizes of approximately 0.1 and 
1 meter, respectively, at the radial location of Jupiter in the minimum mass solar nebula 
model. 



3.2. Initial Setup 

All of our simulations start with uniform gas density (po) and gas pressure (-Pq)- The 
global gas pressure gradient is set to be /3 = —0.04. The initial gas angular velocity profile 
is given by unperturbed Keplerian flow {—{3/2)xQ) and a small reduction due to global gas 
pressure gradient (/3cs/2), Uy = —{3/2)xQ + (3cs/2. Initial disturbances are given to the gas 



- 8 - 



radial velocity with the amplitude \5ux\ = O.OOlCs. Test particles are initially orbiting at 
the Kepler angular velocity. The particles move inward in the initial stage because they lose 
their angular momentum by the headwind that they feel due to the global pressure gradient. 
This is the dust infall problem described in Introduction. Initially 8 particles per gird are 
distributed. Since typical grid we use is (250 - 950) x 100 x 50 (Tabled, totally ~ O(IO^) 
particles are distributed. 

The non-uniformity in the MRI growth is set either by non-uniform resistivity (CASEl) 
or non-uniform vertical (z) component of the magnetic field (CASE2). We describe these 
cases in details in the following. 



3.2.1. CASEl 

In CASEl, gas ionization degree is set to be radially inhomogeneous under the uniform 
magnetic field. The linear analysis for MRI caused b y vertical magnetic field shows th at the 



larger resistivity makes unstable wavelength longer (lJirull996l : ISano &: Miyamalll999l ). The 
upper limit to the wavelength, A^^crit, is scale height for actual disks while in our simulation, 
it is the size of the simulation box in the z direction {Lz). The critical value of resistivity 
with which the modes with wavelength shorter than A^^crit cannot grow is 



1 2 f3Q^ 2k%,,, \ 
?/crit ^ T \ -n - ir^ — H'n, (11) 

where fcz^crit is the wave number corresponding A^^crit- Here Ppiasma = Sc^/f^^ is the plasma 
beta. In our case there is no growing mode if the resistivity makes A^^crit larger than L^. 

In CASEl, we set Ppiasma = 4000 and Bq = (0,0, i?^). Resistivity varies in the radial 
direction such that MRI grows only in a limited zone in the center of the simulation box. 
The resistivity distribution is not changed throughout runs. The radial distribution is shown 
in Figure [2^. In this paper, the zone at the center where resistivity is initially sub-critical 
is referred to as the "unstable" region. The radial width of the stable and unstable regions 
are denoted by and Lu, respectively. 



We have performed only one simulation run for CASEl, which will be called model-?] 
in the later sections. Table [T] lists the parameter values used in the run. 
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3.2.2. CASE2 

The set-up of CASE2 is the same as that in Paper I except for three-dimensional 
simulations. The magnetic resistivity is uniform but the vertical component of magnetic 
field varies in the radial direction. That is, the initial magnetic field is situated as Bq = 
(0, Bq sin 0, Bq cos 6'), where 9 = 9{x) is the angle between the vertical axis and the magnetic 
field. The vertical magnetic component, or 6, determines whether the instability grows or 
not. For large 6, vertical magnetic field is weak and short vertical wavelength modes grow 
faster. Since such modes are dissipated by the ohmic dissipation effectively, the vertical field 
MRI does not develop. The MRI caused by the azimuthal magnetic field is known to grow 



rapid ly with the limit of extremely short vertical wavelength in ideal MHD (IBalbus fc Hawley 



19981) but the such short vert ical wavelength modes tend to be damped resistively. In fact 



Papaloizou &: Terquem (jl997 ) found that the azimuthal field MRI growth ceases for the ion- 



ization degree higher than that stabilizes the vertical field MRI. We also did not find the 
development of the azimuthal field MRI in our simulation as long as the box size in the y 
direction is expanded by a factor up to 3 from the nominal cases at least for about 20 orbits 
in which the gas velocity and dust density vary dramatically. 

The radial distribution of 9 is shown in Figure |2]d. Here, Ppiasma = 400, and the magnetic 
resistivity is ?7 = 0.002H^Q. The larger vertical magnetic component in the center of the 
box makes MRI locally unstable. Outside the unstable zone the value of 6 is smaller than 
the critical one and MRI is stabilized. 

The two-dimensional version of CASE2 was studied in Paper I, in which we found that 
if the spatially averaged magnetic Reynolds number -Rm.ave (Equation ([T])) is small enough, 
the quasi-steady angular velocity profile that is deviated significantly from the Keplerian 
is generated. Since it is expected that this feature is preserved in three dimensional sim- 
ulations and -Rm,ave depends on the radial widths of the unstable region Lu and the stable 
region Lg, we investigate several cases with various Lg while keeping Lu = 1.4H in CASE2: 
Lg = 4.0i/(i?m,ave = 0.096, model-s40), 1.1H{R^^^^^ = 0.37, model-sll), 0.55i/(/?m,ave = 
0.64, model-s055). Note that two-dimensional calculations with the same Lg and were 
performed in Paper I. 

In all of these cases, the friction time is set to be 1.0/Q. To investigate the effects of 
the friction time, we set it to be 0.1/Q in model-tOl, which has Lu = 1.4H and Lg = 4.0H 
(the same as model-s40). The initial settings are summarized in Tabled! 
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3.3. Estimation of Particle Drift Speed 

When the condition i?m,ave < 1 is met, the non-uniform growth of MRI will create a 
new quasi-steady gas angular flow pattern that should affect the dust particle dynamics 
substantially (Paper 1). At the same time, weak remnant turbulent flow may be superposed 
upon the quasi-steady pattern and would affect the dust motions, with its degree being 
dependent on the spatially averaged magnetic Reynolds number. Here we evaluate these 
two effects on the dust accumulation. We will be focusing on the x component of the dust 
velocities Vx- 

If the turbulence is small enough, as we have mentioned in section [21 the velocity 
difference in the y component (angular velocity) between dust and gas is the determining 
factor for the assemblage of dust particles at the boundary between super and sub-Keplerian 
zones. In order to quantitatively analyze the effect of the remnant turbulence, we divide 
the radial component of the particle velocity into two components: one is that due to the 
non-Keplerian gas motion in the quasi-steady state achieved at the MRI saturation and the 
other is that due to turbulence. The equation of motion for the dust velocity fleld at the 
grid point is 

^ = -2n X V + Sfi^xx - — (v - u) . (12) 
dt Tf 

We set V = v' + V]^cp = v' — (3/2)xf2y where v' is a deviation from Kepler flow. Assuming 
dv' /dt = {d\/dt = — (3/2)i7y) in x and y components and eliminating Vy, the x component 
of the particle velocity is given by 

26Uy Ux 

Vx ^ = T -\ ^ 

TfQ + {TfQ)-^ {TfQf+1 

= Vf + vt, (13) 

where 6uy is the gas angular velocity relative to Keplerian velocity, 6uy = Uy — f^ep and 
Ux is the gas radial velocity due to turbulent viscosity. The component Vi represents the 
particle velocity due to non-Keplerian gas motion in the quasi-steady state, because it does 
not vanish even in the non-turbulence case, = 0. Although turbulence also contributes to 
ff , the contribution due to turbulence can be neglected compared with that due to the quasi- 
steady non-Keplerian gas flow, except for strong turbulence case where the quasi-steady flow 
is not achieved. As we will show below, the component Vt can be regarded as characteristic 
amplitude of the velocity dragged by turbulence, although Vt changes with time due to the 
temporal change of Vx and the assumption of dw' /dt = is no more valid in the turbulent 
case. 
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Later in this paper, Vf and Vt will be calculated by equation f lT3|) . using gas velocity ob- 
tained by the MHD simulation, and they will be presented as a function of x after azimuthal 
and vertical averaging. The ratio of the absolute values of the two is 



\Vt\ 
\Vi\ 





Ux\ 


1 


2\6Uy 


TfQ 



(14) 



The turbulence effect is indicated by this ratio. It is stronger, if i) turbulent flow is larger 
than the angular velocity difference 6uy, and/or ii) the dust size (r/) is smaller. In the 
annulus where dust particles are expected to accumulate, 6uy ~ and the ratio diverges. 
However, we will show below that this divergence is restricted only in narrow regions and 
the effect of turbulence does not prevent dust accumulation. 

We have set t he global gas pressure g radient to be /3 = —0.04 so that a direct comparison 



with the result by IJohansen et al.l (120071 ) could be performed. This global pr essure gradient 



however, is less steep than in the Minimum Mass Solar Nebular model ( iHayashil Il98ll ). 
A steeper global pressure gradient {(3 < —0.04) has at least three effects in the current 
context: (1) The baseline of gas velocity becomes lower than the Keplerian velocity and 
it becomes more difficult to create super-Keplerian regions whose outer-edges accumulate 
the dust particles. (2) If a super-Keplerian region appears, accumulation towards its outer- 
edge is faster (t>f is faster) for steeper pressure gradient. This has a positive impact on the 
dust accumulation efficiency. (3) The outer-edge of super-Keplerian region would change 
its location according to the value of (3. The turbulence level at the focal annulus of dust 
accumulation varies accordingly. If the turbulence level is higher, it by itself has a negative 
impact. In CASE2, we have also carried out the runs with (3 = —0.10 and found that the 
trapping of dust particles does occur and its efficiency is slightly higher. 



4. Result 

4.1. Result of CASEl 

In the model-?], the non-uniform excitation of MRI is realized by non-uniform resistivity 
while the magnetic field is set uniform. The result of model-r; nicely illustrates our scenario 
for dust accumulation. 

Figure shows the evolution of MRI. The black lines depict the magnetic field lines 
and the gray scale shows the gas radial velocity. The unstable region lies between two 
white lines {\x/H\ = 0.18). MRI is first excited only in the initially unstable region (see 
the plot at tQ = 19) and significant angular momentum and mass are transported there. 
The MRI turbulence intrudes into the stable region {tQ = 30). Deep inside the stable 
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region, however, is always undisturbed {\x/ H\ > 0.5) because of the rapid dissipation by the 
enhanced resistivity. 

Figures and show the radial profiles of the pressure and the angular velocity of 
the gas, respectively. The quantities have been averaged azimuthally and vertically. The 
sampling times are tVt = 0, 40, and 70. The zone between the two vertical dotted lines is 
the unstable region. The inhomogeneous MRI growth creates the rigid rotation pattern in 
the middle of the simulation box {\x/H\ < 0.15). The pressure distribution is considerably 
modified such that the resultant pressure gradient force balances with the modified Coriolis 
force. The flattened rotation profile cannot sustain the excitation of MRI in the unstable 
zone. Indeed, the turbulence weakens extremely at tQ = 70 in Figure [3^. The profiles of 
pressure or angular velocity in Figures and Hb depict the quasi-steady state set up by the 
non-uniform MRI activity. Most of what we see here is quite similar to the results of the 
two-dimensional simulations described in Paper I even though the initial settings for seed 
magnetic field and resistivity distribution are totally different. 

The rigid rotation causes gas to rotate faster than Keplerian velocity in 0.0 ^ x/H < 0.4 
(FigureHb). This can change the particle migration drastically. FigurelSb shows the temporal 
evolution of the particle density. The color code is set such that the maximum is ten times 
the initial value. After the particles are swept out of the unstable region by the MRI flow, 
they accumulate to the location at x/if ~ 0.4 (note that particles leaving the simulation box 
from the left hand boundary reenter from the right hand boundary after the shearing box 
correction is taken into account). Though not visible in the panels, the particles initially in 
the stable zone are swept likewise towards the same location. The accumulation of particles 
is most clearly shown in Figure Ht in which the radial distribution of the number of particles 
that is averaged azimuthally and vertically and is normalized by the initial value. 

To analyze the particle concentration dynamics in more details, in Figure IHi, we plot 
the maximum (the solid line) and the minimum (the dashed line) values of Vf at a given 
X. The radial velocities of particles are predicted by equation f lT3|) using simulated gas 
velocity after the establishment of the quasi-steady flow {tQ = 70). Since the gas velocity is 
influenced by remnant turbulence, which has variations dependent on y and z, the predicted 
radial velocities have the maximums and the minimums. Figure Hll shows that the turbulent 
fluctuations are sufficiently small compared with the systematic angular velocity change due 
to the modulation of the pressure profile. In Figure Htl, the max/min of radial amplitude 
of turbulent velocity, Vt, is plotted. The amount of ft is smaller than those of ff except 
in the region l^fl ^ 0. Thus, Vf represents the typical radial velocity of particles that is 
induced by head/tail angular wind in the absence of turbulence. Also plotted by dots are 
the actual particle radial velocities in the simulation result, Vx, which are concentrated near 
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the super/sub-Keplerian boundary. The data are taken at tQ = 70, but it stays essentially 
the same after the MRI saturation at tfi ^ 40. In the radial regions where the maximum 
ff is negative, gas rotation is sub-Keplerian for all y and z. In such radial locations, all 
the particles lose their angular momentum by the drag due to the slower angular wind 
and migrate inward. Conversely, in the regions where the minimum Vf is positive, all the 
particles migrate outward. In the middle of the stable zone {\x/H\ > 0.5), the gas flow is 
hardly changed. The value of V{ ~ — 0.02cs corresponds to the infall speed due to global 
pressure gradient under the initial condition. 

The particles are concentrated to the zone where both min[ff] < and max[ff] > 
are satisfied, which is situated around x/if ~ 0.4. The zone satisfying the two inequalities 
is rather narrow in x, and the difference between the maximum and the minimum Vf in 
this zone is small, resulting in high concentration of the dust particles at the outer-edge of 
the super-Keplerian zone as shown in Figure Ht. Because \vt/{vf + tit)] ^ 1 in most of the 
regions except for the dust concentration zone, and the width of the dust concentration zone 
is small, the effects of turbulence do not significantly expand the dust concentration zone. 

In Figure we plot the time variation of the maximum number of particles per grid. 
The more or less monotonic increase leads to the peak density of more than 1,000 times the 
initial density for tfi > 50. We will show that the same clump grows in density while keeping 
its identity, which is needed for subsequent gravitational instability. 

We pick up a particular time tc during the period when the number of dust particles 
in the densest grid in the whole region increase monotonically or is saturated. Then, we 
search for the cell which has the largest number of particles in it. This cell is identified 
as the center of the most prominent clump that is composed of a number of cells. The ID 
numbers of the Nc particles in this cell of the highest density at t = tc are recorded and 
their motions are traced in time (also backward in time if necessary). To determine a new 
position of the center of the same clump at different times, we inspect the cells within 5 
grid distance {0.05H) from the center of mass of the traced particles. The cell having the 
largest number of particles among them is defined to be the new position of the center of the 
clump. Table [2] shows t^. and for the runs described in this paper. The time evolution of 
the number of particles in the cell at the center of the traced clump is plotted by the dotted 
lines in Figure [5^. It agrees with the maximum number of particles over the whole region, 
indicating that the same clump keeps its identity and stays to be the most prominent clump. 

Another way to confirm that the clump steadily holds together is to inspect the velocity 
dispersion of particles. Figure presents the velocity dispersion of particles as a function of 
distance from the center of the clump ait = tc- All the compositional velocity dispersions are 
extremely low (< 0.5m/s), indicating that the particles are not to be diffused significantly. 
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When the self-gravity is taken into account, the high mass density and the low velocity 
dispersion should be the preferable condition for the subsequent planetesimal formation via 
gravitational instability. 

When the global pressure gradient is steeper, the effects on dust dynamics described at 
the final part of Section 3 must be considered. For the case shown here, we make discussion by 
inspecting Figure HJi and equation (fT3l) . From the maximum amplitude of min[ff] change (~ 
O.OScs), we can expect that the value of min[ff] remains positive (vertically and azimuthally 
uniform super-Keplerian region survives) if (3 > —0.16. That is, dust particles are prevented 
from infall as long as /3 > —0.16. This condition would be marginally satisfied at r < 5AU 
in the MMSN model. Though the dust concentration that occurs at f f ^ would be shifted 
inward with the more downward shift of the Vf (larger negative the turbulence effect on 
dust concentration would not change because the turbulence level Vt is evenly low. Therefore, 
with a steeper global pressure within the range of (3 > —0.16, the dust particles would drift 
faster towards and accumulate more efficiently at the outer-edge of super-Keplerian region. 

4.2. Result of CASE2 

In the models of CASE2, the non-uniform excitation of MRI is realized by non-uniform 
distribution of the magnetic field while resistivity is set uniform. In the presence of the 
substantial dissipation, the variation of the vertical component of the magnetic field switches 
the MRI condition from stable to unstable. In CASE2, the magnetic field angle is varied 
along the x (radial) direction such that MRI unstable region is situated only in the middle 
of the simulation box. The overall MHD dynamics varies according to the spatially averaged 
magnetic Reynolds number (-Rm,avc) which is controlled by the width of the unstable/stable 
regions in the simulation box, as shown in Paper I. Here, the rotation pattern of the magnetic 
field and the width of the unstable region (Lu = lAH) is fixed. 

4.2.1. model-s40 

In the model-s40, the radial width of initial stable region Lg is 4.0H. The relatively 
wide stable region results in -Rm,ave = 0.096(<C 1), which ensures effective enough dissipation 
to realize the quasi-steady state. Figure |6^ presents the evolution of MRI in model-s40. 
The gray scale and lines are radial velocity of gas and magnetic field lines, respectively. 
The boundaries between the stable and unstable regions are depicted by the two white lines 
{\x/H\ ^ 0.7). It is found that the results of the present three-dimensional run are similar 
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to the corresponding two-dimensional results described in Paper I. The instability is excited 
only inside the unstable region. A part of the stable area is disturbed as well, because the 
resistivity is not large enough to dissipate the fluctuations immediately {tQ = 35). They 
are, however, eventually dissipated (tQ = 40 and 80). 

Figure [7| presents radial distributions of pressure and angular velocity of gas, and the 
vertical component of the magnetic field. All the quantities obtained at tQ = 0, 40 and 
70, are averaged both azimuthally and vertically. The inhomogeneous MRI growth creates 
the near-rigid rotation pattern in the middle {\x/H\ < 0.6) and the modified Coriolis force 
balances with the newly created pressure gradient force. The small vertical magnetic field 
together with the fiat angular velocity profile suppresses the growth of MRI at later times 
inside the unstable region in which MRI was initially excited. 

While these results of model-s40 are similar to those of model- (CASEl), three differ- 
ences are pointed out: i) The pressure and angular velocity are distributed less smoothly in 
the stable region in Figures [7^ and[7)D, which is attributed to deeper penetration of the MRI 
activity, ii) super-Keplerian region exists at x/iJ ^ —2.5, in addition to the middle region 
of 0.0 ^ x/H < 2.0, and iii) the weak instability persists at least until tQ = 90 in the stable 
region near the unstable region (Figure [6^). The super-Keplerian region at x/H ^ —2.5 
arises, because larger stable region makes mass and angular momentum transport less ef- 
ficient around the middle of the stable region [x/H ~ ±4.5). As a result, the pressure 
gradient becomes positive at a;/-ff ~ —2.5. In model- r/, the large magnetic resistivity in the 
stable region suppresses the instability even in the region of steep radial gradient of angular 
velocity. But, since in the models of CASE2, constant magnetic resistivity is assumed, the 
magnetic field is not completely dissipated in the stable region with the steep radial gradient. 

Figure presents the temporal evolution of the particle density. The evolution of the 
particle density pattern is similar to CASEl. However, unlike CASEl, not all the particles 
are assembled to the focal annulus. Figure [8^ shows the azimuthally and vertically averaged 
value of the number of particles as a function of x. As expected, particles are condensed at 
x/if ~ 2.0 (the outer-edge of the super-Keplerian region). Figure [8^ also shows some other 
peaks in the stable region, which reflects the characteristics i) and ii) of the flow pattern 
described above. 

Figure [8]d shows radial velocities of particles predicted by equation (fT3|) using simulated 
gas velocity after the establishment of the quasi-steady flow {tQ = 70). The dots concen- 
trated at x/H ~ 2.0 in the figure are actual velocities of all the simulated particles. The 
actual data of the particle velocities are well reproduced by equation (IT^ , since most of the 
data points are located between the predicted minimum and maximum values. We further 
decompose Vd into Vf and Vt according to equation f|T3l) . Figure (Ht presents the radial distri- 
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bution of the two components. The maximum and minimum values for each component are 
plotted. It is shown that the effect of the remnant turbulence is small and overall features of 
particle radial migration is represented by V{. Due to the non-smooth pressure and velocity 
distributions, min[fd] — min[t'f] > at —3.5 ^ x/H < —2.5 and x/H ^ —2.1, in addition 
to 0.0 < x/H < 2.0. It means that all particles have positive radial velocity there and 
accumulate near the outer edges of these regions. On the other hand, max[wd] — max[ff] > 
is observed in some areas, for example, around x/iJ ~ 3.0. A small number of particles are 
stalled there. 

The peaks of particle density at a;/if < would not grow because there are few particles 
which can accumulate from outside (Figure [8^). The stagnant areas cannot stall all the 
particles which migrate inward due to min[wd] < (Figure [Hb). Therefore the peak at 
x/if ~ 2.0 is the most prominent concentration zone. 

Figure [9^ shows the time variation of the number of particles in the densest grid in 
the whole region and at the center of the traced clump. The agreement between the two 
indicates that the same clump is growing. Figure [Hb is the velocity dispersion of particles 
within the clump at t = tc(= 55/Q). The velocity dispersion is small but is larger than in 
model-?]. While the turbulence does not affect the location of the dust concentration, it does 
influence the degree of particle accumulation. 

When the azimuthal box size is loi ig, shear velocity be c omes supersonic near the edge 



and it may cause artificial density dip (jjohnson et al.l l2008l : iJohansen et al.l 120091 ) unless a 



high-order scheme is applied for Keplerian advection term. To confirm that the density dip 
found in our simulations is not caused by such a numerical error, we also carried out a run 
in which the unstable region is put in the side areas and the stable region is at the center. 
The other parameters are the same as model-s40. We found that the density dip is formed 
at the side area and no artificial density dip is created in the central area, which means that 
our numerical scheme using CIP method with the staggard mesh and the short timestep 
{dt = 10~^; Courant condition is set to be < 0.5) is sufficiently high-order scheme. 

Let us discuss the effects of a steeper global pressure gradient. Following the same 
argument as in the final part of section 4.1, from the value of Vf or 6uy, a super-Keplerian 
region would remain unless the global prssure gradient is as anomalously larg e as /3 < —0.64. 



To make a more quantitative argument, as was done in IJohansen et al.l (120071 ). we also carried 
out a run with the pressure gradient of /3 = —0.10 (model-s40b). The other parameters are 
the same as model-s40. Figure [TUl shows the close-up of the distribution of max, min[ff] and 
max, min[ft]. With the higher negative value of (3, gas rotates more slowly and the stronger 
headwind deprives the dust particles of angular momentum more efficiently to accelerate 
their inward migration in the sub-Keplerian region. Since ft does not depend on /3, the value 
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of |ff I exceeds |ft| by a larger degree in outer regions. Thus greater amounts of dust particles 
move to the dust concentration area without becoming locally stagnant due to non-smooth 
gas pressure and gas velocity fluctuations in the sub-Keplerian region. On the other hand, 
in the dust concentrated area, the dust particles are affected by turbulence more severely. 
With more downward shift of the Uf value in the whole region, the dust concentration that 
occurs at Uf ~ is shifted inward to x/if ~ 1.8 from x/H ~ 2.0. The magnetic turbulence is 
stronger at this new location and it increases the radial width {5x) of the dust accumulated 
region (denoted by max[fd] > and min[fd] < 0: the gray hatched areas in the figure from 
5x/H ~ 0.5 to 5x/H ~ 0.9 and doubles the velocity dispersion of particles there. To see 
the overall consequences of these both positive and negative effects, the number of particles 
in the most dense grid cell is plotted in Figure [9^ together with the result of model-s40. 
Despite the different dust radial velocity and gas turbulence, the density increase rate and 
the peak value are only slightly different from each other, although the new run could show 
higher dust density on longer timescale. 



4.2.2. model-to 1 

The effects of turbulence are larger for particles with shorter tj. We have run a case with 
TfVL = 0.1 (model-tOl). Except for the friction time, all the settings are the same as model- 
s40. Since there is no kick-back from the dust particles to the MHD fluid, MHD results are 
exactly the same (Figure [6^ and Figure [7]). The difference in the dust concentration pattern 
is caused by the enhanced gas drag. 

Figure presents the evolution of the particle density. The difference from the previous 
case that becomes visible at the later time is that the smaller dust particles are scattered 
over a broader area. Figure [TTh shows the radial proflle of the dust number density averaged 
azimuthally and vertically. Since the MHD flow patterns are exactly the same, the prominent 
density peaks are located at the same position as in model-s40. However, the peak values are 
lower and the radial widths are wider. Furthermore, more particles stay outside the peaks 
and are scattered in the stable region. 

These differences are caused by the enhanced turbulent diffusion of dust particle motions 
through enhanced gas drag. Indeed, Figure [TTb shows that the turbulent scattering ft 
dominates over the steady migration Vf in most areas of the simulation box. Figure [TTb 
shows the estimated radial velocity of particles Vd- The radial velocity of the simulated 
particles, Vx, are also shown by dots. Although dust particles are swept out of the unstable 
region, in most parts of the unstable region, max[fd] > and min[fd] < are satisfled, 
so the dust particles are broadly distributed in the stable region with lower concentration. 
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Nevertheless, identity of the clump with the highest dust concentration is still maintained 
(Figure [T2b). Indeed the velocity dispersion around center of the traced clump shown is not 
increased significantly (Figure [T^b)- 

We also run a high global pressure gradient model (/3 = —0.10). Though the inward 
migration of dust particles becomes faster on average, the maximum density is unchanged 
because the remnant turbulence that scatters the dust particles in the dust concentrated 
region is the dominant factor for this tjQ = 0.1 case. 



4.2.3. model-s055 



Recent works (e.g., IJohansen fc YoudinI 120071 ) proposed that turbulence itself concen 



trates dust particles in turbulent eddies. In their models, lifetime of the eddies has to be 
longer than the timescale to form dense enough particle clumps for formation of planetes- 
imals. On the other hand, our model proposes another path to accumulate dust particles 
that turbulence plays a role in transformation of the nearly Keplerian gas flow into the quasi- 
steady flow with a local rigid rotation region. Because the dust particles accumulate near 
the outer edge of super-Keplerian parts produced by the rigid rotation and the flow pattern 
is quasi-steady, the timescale problem does not exist in our model. Actually, in the case of 
model-s40 and model-tOl with -Rm.ave — 0.096, the same clump with the highest density is 
maintained until the end of simulations {tQ ~ 80). 

However, in our model, as -Rm,ave increases to unity, stronger residual turbulence may 
destroy the clumps, although the clumps are repeatedly created, which may inhibit plan- 
etesimal formation. Here, we examine the cases with larger -Rm,ave (but still -Rm,ave < 1) to 
find the critical value of -Rm.ave for persistent clumping. 

In model-s055, the width of the stable region Lg is one eighth of model-s40. Accordingly, 
-Rm,ave = 0.64 Compared to 0.096 of model-s40. Figure [T5k shows that the MRI turbulence 
intrudes into the stable region and covers the whole region without being dissipated. By 
tQ = 80, the instability is weakened. Although a clear rigid rotation is not formed inside 
the unstable region (Figure [T4k ) because of the stronger turbulence that remains especially 
in the stable region, the velocity gradient becomes small inside the unstable region and the 
super-Keplerian zone is formed. As a result, the dust particles are swept out of the unstable 
region and accumulate near the outer stable/unstable boundary (Figure IT^). 

However, because the unstable region occupies larger fraction of the simulation box, the 
resultant magnetic field diffused out of the unstable region is larger than in the previous 
models. These lead to stronger remnant turbulence. Figure fTEk shows that the number of 
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particles in the densest grid of the whole region does not agree with that at the center of the 
traced clump. It means that the clumping is only tentative, although clumps are repeatedly 
created. The diffusive nature of the clumps in model-s055 is depicted by the larger velocity 
dispersion shown in Figure fTEb . 

Thin hues in Figure [T5k shows the clump identity test for model-sll (-Rm.ave = 0.37). 
It shows that the clumps are not persistent as well. Therefore, we conclude that persistent 
clumps of dust particles are formed for -Rm,ave ^0.1. 

If the global pressure gradient is steeper, the radial velocity of dust particles is faster 
and its effect may seem to compete better with that of gas turbulent motion. However, from 
Figure [T3Ji, the minimum value of V{ becomes negative in the entire regions if /3 < —0.24 and 
no particles can be accumulated at the boundary of super/sub-Keplerian regions. There is 
not much parameter space where more dust concentration can be expected. 

5. Conclusion and Discussion 

We have investigated the dust concentration in protoplanetary gas disks that are in 
the quasi-steady state created by inhomogeneous magnetorotational instability. We set the 
inhomogeneous instability with the initial radial non-uniformity of either ionization degree or 
vertical component of the magnetic field. The result of our local three-dimensional resistive 
MHD simulations with Lagrangian particles clearly shows that the meter or decimeter sized 
dust particles are concentrated strongly and steadily in either setting except for the cases 
with -Rm,ave > 0.1 (where -Rm.ave IS the averaged magnetic Reynolds number averaged over 
the simulation box), which can allow planetesimal formation via gravitational instability. 

The MRI growth in our model generates the locally confined super- Keplerian region and 
this flow is quasi-stable due to the support of gas pressure gradient which is also formed by 
the MRI. Because unperturbed gas flow is sub-Keplerian due to global pressure gradient, 
the particles suffering gas drag are concentrated in the Keplerian domain at the outer-edge 
of the super-Keplerian area. Since the flow is in a quasi-steady state, particles are supphed 
from outer regions by gas drag migration and the particle density increases significantly. 

The process of particle concentration and the increasing rate of dust density depend 
on the particle size, the initial radial widths of unstable/stable regions, and the initial non- 
uniformity setting as follows: 

The initial non-uniform setting, magnetic field or resistivity: The non-uniform 
resistivity model produces more stable and clean state of local rigid rotation than the 
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non-uniform magnetic field model, because the high resistivity dissipates the magnetic 
perturbations more rapidly in the stable regions. Indeed, the velocity dispersion of 
particles in clumps is smaller in the non-uniform resistivity model. In the non-uniform 
magnetic field cases, the magnetic reconnection and the remnant weak instability create 
less smooth gas velocity field, so that particles tend to assemble at multiple radial 
positions. On the other hand, in the non-uniform resistivity case, all particles are 
concentrated at a specific radial position, and the peak values of the particle density 
in clumps is higher, although the particle density also depends on the size of numerical 
box. 

The dust size, meter or decimeter: The meter-size particles are less coupled with gas 
motion than decimeter-size ones. Since their radial migration is faster, the meter- 
sized particles are more likely to be condensed in the radially narrow area without 
being scattered by the remnant gas turbulence. The peak dust density becomes > 
1000 times larger than the initial value. Even in the decimeter-size particle cases, 
the velocity dispersion of particles in the clump excited by the remnant turbulence is 
sufficiently smaller than radial drift velocity due to the quasi-steady gas flow and dust 
concentration is still observed in the simulations, although the enhanced dust density 
is at most ~ 100. 

The degree of remnant turbulence: When the stable region is initially smaller in the 
non- uniform magnetic field model (CASE2), that is, -Rm.ave is larger, stronger turbu- 
lence tends to remain. For -Rm,avc ^ 1, however, the super-Keplerian area continues 
to exist and clumps with density ~ 100 times larger than initial values are robustly 
created. Furthermore, for -Rm,ave ^0.1, the clumps are persistent enough for following 
gravitational instability, while they are tentative and repeatedly produced in the cases 
of0.1<i?„<l. 

We should address the following points in subsequent papers: 

consistent magnetic resistivity: We need to calculate the value of magnetic resistivity 
consistently as a time dependent value, because it is sensitive to both gas and dust den- 
sities. Around strong dust concentration areas, gas and dust density become larger and 
magnetic resistivity becomes higher, which may calm down the remnant turbulence. 

Rossby wave instability (RWI): RWI can be caused from a pressure bump on the radial- 
azimuthal plane. If we simulate in an enlarged box, vortex may grow, although in the 
non-uniform magnetic field cases, the growth of the instability may be reduced by 
the azimuthal magnetic field. If RWI o ccurs, it may serve as a mechanism o f dust 



concentration in the azimuthal direction (jinaba fc Bared l2006l : iLyra et al.ll2008l ). 
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the vertical structure of a protoplanetary disk: In the upper disk regions, ionization 

degree may be so high t hat the MRI to occur (jFleming fc Stone 2003 : Fromang fc Papaloizou 



20061 : lOishi et al.ll2007t ) . The turbulence there may affect the small dust dynamics near 



the midplane. 

feedback from dust particles onto gas: The feedback may be the most important. A 
great concentration in our results here will inevitably change gas velocity field and it 
could affect the condition of dust concentration (density enhancement, velocity disper- 
sion and so on) and planetesimal formation. We have already started simulations with 
this effect. 

collisional destruction: While the velocity dispersions of dust particles in clumps are 
extremely low in some of our results, the dust particles can be collisionally disrupted 
into small fragments that are coupled to gas motion and the fragments may be dispersed 
by even weak turbulence. For clumps to gravitationally collapse, its timescale needs 
to be smaller than mean collision time. On the other hand, if the dust grains are 
also produced by the collisions, they might affect the gas ionization state through 
recombination of electrons onto their surface. We need to address these feedback 
effects by collisional disruption as well. 
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Fig. 1. — Results of the two-dimensional version of model-s40 {dBz{t = 0)/dx ^ 0,Ls/H = 
4.0) described in Paper I. The corresponding result of the three-dimensional simulation is 
shown in Fig. 7. Time evolution of vertically averaged values of (a) pressure P and (b) 
gas angular velocity Uy. The thin solid, dotted, dashed and bold lines are the snapshots at 
tfl = 0, 27, 40 and 70, respectively. The two vertical dotted-lines are the boundaries between 
the unstable and stable regions. MRI is initially excited only between the two dotted lines. 
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Fig. 2. — Two versions of initial setting to realize the non-uniform MRI growth in the 
simulation box. (a) The non-uniform radial distribution of resistivity used in CASEl (model- 
T]), where the dashed line represents the critical value for MRI development predicted by the 
linear theory, (b) The non-uniform radial distribution of the magnetic field angle 6 in CASE2 
(model-s40), where the dashed line shows the critical angle for MRI predicted by the linear 
theory (for detail, see Paper I). 
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Fig. 3. — The results of model-77 {drj/dx 7^ 0, TfQ = 1.0). (a)Time evolution of the magnetic 
field (black lines) and the angular velocity Uy (contours) of gas. (b) Particle density normal- 
ized by the initial value. The unstable region is between the two white lines. Non-uniform 
growth of MRI and its relaxation seen in Panel a lead to the dust concentration depicted in 
Panel b. 
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Fig. 4. — Radial dependences of (a) pressure P, (b) gas angular velocity Uy, (c) the number 
of particles normalized by the initial value, and (d) the maximum and minimum values of 
ff and ft and the actual particle radial velocity v^^i, in model-r] (dr^/dx 7^ 0,TfQ = 1.0). 
They are averaged both azimuthally and vertically. In Panels a, b and c, the dotted, dashed, 
and solid lines represent the results at tQ = 0,40 and 70, respectively. Since V{ is predicted 
by equation f|T3|) using simulated gas velocity after the establishment of the quasi-steady 
flow {tfl = 70) and the simulated gas velocity has a component of turbulence, it has the 
maximum (the solid line) and minimum (the dashed line) values as shown in Panel d. In 
Panel d, the maximum/minimum values of Vt and the actual particle velocity are also plotted 
by dash-dotted/dotted lines and dots. 
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Fig. 5. — Time evolution of dust concentration in model-77 (d?7/dx 7^ 0,TfQ = 1.0). (a) The 
solid line represents the number of particles in the cell having the highest density in the 
whole region, which is normalized by the initial value, while the dashed line is the number of 
particles in the cell at the center of the traced clump, (b) Velocity dispersion in the traced 
clump at tQ = t^^l (=55.0 in this case). The solid, dashed, dash-dotted and dotted lines 
show total, radial, azimuthal and vertical components, respectively. The sound velocity is 
estimated as 500 m/s with temperature T ~ 70K at 5AU. 
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Fig. 6. — (a)(b) The same as Figure [3] but for model-s40 {dBz(t = 0)/dx ^ 0,Ls/H = 
4.0, TfQ = 1.0). (c) Particle density distribution for model-tOl {rfQ = 0.1). 
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Fig. 7. — Time evolution of azimuthally and vertically averaged values of (a) pressure P, (b) 
gas angular velocity Uy, and (c) vertical magnetic component B^, for model-s40 and model- 
tOl {Ls/ H = 4.0). The dotted, dashed and bold lines represent the snapshots at tfl = 0,40 
and 70, respectively. 
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Fig. 8. — Dust concentration in model-s40 {Ls/H = 4.0, TfQ = 1.0). (a) Time evolution 
of the number of particles in each radial grid scaled by the initial value, where the dotted, 
dashed and bold lines represent the snapshots at tfl = 0,40 and 70, respectively, (b) Radial 
velocity of particle (dots) and the estimated radial velocity that are azimuthally and 
vertically averaged, where the solid and dashed lines express the maximum and minimum 
values in the averaging, respectively, (c) ff (max: bold solid, min: bold dashed) and max/min 
of Vt (thin lines). The results in Panels b and c are obtained at tQ = 70. 
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Fig. 9.— The same plots as Figure [5]but for model-s40 {LJH = 4.0, r/fi = 1.0, /3 = -0.04). 
The dotted hne in Panel a is the number of particles in the densest cell in model-s40b 
{/3 = -0.10). 
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Fig. 10. — (a) The closeup of Figure E}:; (model-s40 in which /3 = —0.04). (b) The same plots 
with Panel a but for model-s40b in which /3 = —0.10. 




Fig. 11. — The same plots as Figure |8] but for model-tOl {L^/H = 4.0, TfQ = 0.1), where 
turbulence effects are enhanced. 
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Fig. 12. — The same plots as Figure |5]but for model-tOl {L^/H = 4.0, TfQ = 0.1). 
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Fig. 13. — The same plots as Figure |3] but for model-s055 {Ls/H 



= 0.55, TfQ = 1.0). 
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Fig. 14. — (a) Time evolution of both azimuthally and vertically averaged angular velocity 
of gas Uy in model-s055 {L^/H = 0.55, TfQ = 1.0), where the dotted, dashed and bold 
lines represent the snapshots at tfl = 0,40 and 70, respectively. (b)(c)(d) The same plots 
as Figure [Hk, b, and c but for model-s055. Dust particles are scattered due to elevated 
turbulence. 
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Fig. 15. — The same plots as Figure |5] but for model-s055 {Ls/H = 0.55, r/^7 = 1.0). The 
thin dash-dotted and thin dotted hues in Panel a show the results of model-sOll {L^/ H = 

1.1, Tfil = 1.0). 
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model- 


L^x LyX {/H) 


Nr,X NyX 


non-uniform 
parameter 


Lu/H 










V 


1.4 X 1.0 X 0.14 


280 X 200 X 28 


V 


0.28 


0.56 


*** 


1.0 


1.3 


s40 


9.5 X 1.0 X 0.50 


950 X 100 X 50 


B 


1.4 


4.0 


0.096 


1.0 


3.8 


sll 


3.5 X 1.0 X 0.50 


350 X 100 X 50 


B 


1.4 


1.1 


0.37 


1.0 


1.4 


s055 


2.5 X 1.0 X 0.50 


250 X 100 X 50 


B 


1.4 


0.55 


0.64 


1.0 


1.0 


toi 


9.5 X 1.0 X 0.50 


950 X 100 X 50 


B 


1.4 


4.0 


0.096 


0.1 


3.8 



Table 1: RUN PARAMETERS (1): Name of run. (2): Size of simulation box. (3): Grid 
resolution. (4): Parameter set non-uniformly. (5): Radial width of unstable region. (6): 
Radial width of stable region. (7): Friction time. (8): Number of particles. 
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model- 






Nc 




V 


55.0 


34.0 


71096 


0.013(0.0023, 0.012, 0.0036) 


s40 


75.0 


37.4 


13823 


0.16(0.041, 0.14, 0.054) 


sll 


75.0 


33.6 


2285 


0.77(0.67, 0.24, 0.13) 


s055 


90.0 


43.2 


392 


0.65(0.57, 0.24, 0.13) 


toi 


75.0 


11.2 


292 


0.18(0.086. 0.13. 0.069) 



Tabic 2: TRACED PARTICLES (1): Name of run. (2):Timc to define which particles are 
traced. (3): Period for tracing. (4): Number of traced particles. (5): Position dispersion of 
Nc particles averaged for 5ttrace^ and its composition. 



